if !d.name eq 'PS' then begin
    device,xsize=18,ysize=12,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end
default,png,0
window,1,xsize=480,ysize=640
!x.margin=[10,3]
!p.charsize=1.5
; field to be ploted
@eu_setup
field=bz
;;;;;;;;;;;;;;;;;;;;;
default,ibeg,0
nt=(size(field))[3]
print,nt
nl=64
year = 365.*24.*60.*60.
time = (nxaver*dt*indgen(nt))/year
good=where(time gt 50)
rg=where(r gt 0.74) 
bmax=(max(abs(field[*,rg,good])))/0.33
bmin=-bmax
lev=grange(-0.1,0.1,nl)
year = 365.*24.*60.*60.
fo = "(F5.2)"
ii=0

cpos=[0.23,0.08,0.99,0.9]
bpos=[0.37,0.33,0.39,0.655555]
for i=ibeg, nt-1 do begin
 time = nxaver*dt*double(ii)/year
 contour,smooth(transpose(reform(field[*,*,i])),2),x[*,*],y[*,*],$
         /fi,nl=64,/iso,lev=lev,pos=cpos,title='!8B!D!7u!N!8 / B!D!6eq!6!N, !8t=!6'+string(time,FORMAT=fo)+' !6yr'
 colorbar,range=[min(lev),max(lev)],/vertical,xaxis=2,$
         ytickformat='(F6.2)',yticks=4,pos=bpos,charsize=1.5,$
         ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)]
 plots,circle(0,0,0.97),thick=3
 plots,circle(0,0,0.61),thick=3
 plots,circle(0,0,0.718),thick=2,li=2
 oplot,[0,0],[-0.97,-0.61],thick=2,li=0
 oplot,[0,0],[0.61,0.97],thick=3,li=0
 print,i
    if (png EQ 1) then begin 
     istr2 = strtrim(string(ii,'(I20.4)'),2) ;(only up to 9999 frames)
     print,"writing png",'./img_'+istr2+'.png'
     image = tvrd(true=1)
     imgname = './img_'+istr2+'.png'
     write_png, imgname, image, red, green, blue
    endif
 wait,0.1
 ii ++
endfor

end

